{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from astropy.io import fits\n",
    "from astropy.modeling import models, fitting\n",
    "import photutils\n",
    "#import lmfit\n",
    "#from lmfit.lineshapes import gaussian2d\n",
    "#from lmfit import Parameters, Model\n",
    "import numpy as np\n",
    "import math\n",
    "from astropy import wcs\n",
    "from astropy.coordinates import SkyCoord\n",
    "from aplpy import FITSFigure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this script is to identify the shift in WCS positions between SOFIA Flights 684 and 685. \n",
    "\n",
    "1- Load the level 2 data for AOR 92 in flights 684 and 685.\n",
    "\n",
    "2- Fit 2D Gaussian profile to bright source Northeast of the map to find the average WCS coordinate of the center, and then calculate the shift between the two flights.\n",
    "\n",
    "3- Hack the CRVAL1/2 keywords in the header of each extension of the level 2 data in flight 684 using the measured shift."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F0684_HA_POL_08018692_HAWEHWPE_SMP_013-024.fits\n",
      "F0684_HA_POL_08018692_HAWEHWPE_SMP_049-060.fits\n",
      "F0684_HA_POL_08018692_HAWEHWPE_SMP_061-072.fits\n",
      "F0684_HA_POL_08018692_HAWEHWPE_SMP_013-024.fits\n"
     ]
    }
   ],
   "source": [
    "# Reading all the files using AOR 92 for flight 684\n",
    "RefDir684 = 'Reference/684_92/'\n",
    "\n",
    "RefFiles684 = []\n",
    "for entry in os.scandir(RefDir684):\n",
    "    if entry.is_file() and entry.name.endswith('.fits'):\n",
    "        print(entry.name)\n",
    "        RefFiles684.append(entry.name)\n",
    "\n",
    "RefSize684 = len(RefFiles684)\n",
    "\n",
    "print(str(RefFiles684[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F0685_HA_POL_08018692_HAWEHWPE_SMP_001-012.fits\n",
      "F0685_HA_POL_08018692_HAWEHWPE_SMP_013-024.fits\n",
      "F0685_HA_POL_08018692_HAWEHWPE_SMP_025-036.fits\n",
      "F0685_HA_POL_08018692_HAWEHWPE_SMP_037-048.fits\n",
      "F0685_HA_POL_08018692_HAWEHWPE_SMP_001-012.fits\n"
     ]
    }
   ],
   "source": [
    "# Reading all the files using AOR 92 for flight 685\n",
    "RefDir685 = 'Reference/685_92/'\n",
    "\n",
    "RefFiles685 = []\n",
    "for entry in os.scandir(RefDir685):\n",
    "    if entry.is_file() and entry.name.endswith('.fits'):\n",
    "        print(entry.name)\n",
    "        RefFiles685.append(entry.name)\n",
    "\n",
    "RefSize685 = len(RefFiles685)\n",
    "\n",
    "print(str(RefFiles685[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating the box for the analysis\n",
    "BoxSizeX = 0.02 # degrees\n",
    "BoxSizeY = 0.02 # degrees\n",
    "CenterX_684 = 277.7313838 # degrees\n",
    "CenterY_684 = -10.5003242 # degrees\n",
    "CenterX_685 = 277.6822875 # degrees\n",
    "CenterY_685 = -10.3693331 # degrees\n",
    "Angle = 0.0 # degrees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def crop_map(FitsFile, BoxSizeX, BoxSizeY, CenterX, CenterY):\n",
    "    # Reading size of the data container\n",
    "    SizeFile = np.shape(FitsFile[0].data)\n",
    "    # Creating a copy of the container that will be modified\n",
    "    CropFile = FitsFile[0].data.copy()\n",
    "    # Reading the WCS information for the map\n",
    "    WCSinfo = wcs.WCS(FitsFile[0].header)\n",
    "    # Cropping the map to keep only the requested region\n",
    "    for i in range (0,SizeFile[0]): # Dec\n",
    "        for j in range (0,SizeFile[1]): # RA\n",
    "                    # Current pixel position\n",
    "                    PixCoord = SkyCoord.from_pixel(j,i, WCSinfo)\n",
    "                    CurrentRA = PixCoord.ra.deg\n",
    "                    CurrentDec = PixCoord.dec.deg\n",
    "                    # Distance to central pixel\n",
    "                    delta_ra = (CurrentRA - CenterX) * math.cos(CurrentDec*math.pi/180.0)\n",
    "                    delta_dec = CurrentDec - CenterY\n",
    "                    # Coordinate conversion\n",
    "                    x1 = delta_ra * -1.0\n",
    "                    y1 = delta_dec\n",
    "                    ang = Angle*math.pi/180.0\n",
    "                    delta_x2 = abs(x1*math.cos(ang)+y1*math.sin(ang))\n",
    "                    delta_y2 = abs(-x1*math.sin(ang)+y1*math.cos(ang))\n",
    "                    # Is the pixel in the mask?\n",
    "                    if delta_x2 < 0.5*BoxSizeX and delta_y2 < 0.5*BoxSizeY:\n",
    "                        nothing = 0\n",
    "                    else:\n",
    "                        CropFile[i,j] = np.nan\n",
    "    # Creating the new fits file\n",
    "    CropFile_FITS = fits.PrimaryHDU(data=CropFile, header=FitsFile[0].header)\n",
    "    # Returning new fits file\n",
    "    return CropFile_FITS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating list of fits files\n",
    "ListFiles684 = []\n",
    "ListCrops684 = []\n",
    "# Cropping the map for each file\n",
    "for k in range (0, RefSize684):\n",
    "    # Opening FITS file\n",
    "    File684 = fits.open(RefDir684+RefFiles684[k])\n",
    "    # Creating a copy of the file\n",
    "    CopyFile = File684.copy()\n",
    "    ListFiles684.append(CopyFile)\n",
    "    # Creating cropped map\n",
    "    Crop684 = crop_map(File684, BoxSizeX, BoxSizeY, CenterX_684, CenterY_684)\n",
    "    # Creating list of FITS files\n",
    "    ListFiles684.append(File684)\n",
    "    ListCrops684.append(Crop684)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating list of fits files\n",
    "ListFiles685 = []\n",
    "ListCrops685 = []\n",
    "# Cropping the map for each file\n",
    "for k in range (0, RefSize685):\n",
    "    # Opening FITS file\n",
    "    File685 = fits.open(RefDir685+RefFiles685[k])\n",
    "    # Creating a copy of the file\n",
    "    CopyFile = File685.copy()\n",
    "    ListFiles685.append(CopyFile)\n",
    "    # Creating cropped map\n",
    "    Crop685 = crop_map(File685, BoxSizeX, BoxSizeY, CenterX_685, CenterY_685)\n",
    "    # Creating list of FITS files\n",
    "    ListFiles685.append(File685)\n",
    "    ListCrops685.append(Crop685)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Input data contains non-finite values (e.g., NaN or inf) that were automatically masked. [photutils.centroids.gaussian]\n",
      "WARNING: Input data contains non-finite values (e.g., NaN or inf) that were automatically masked. [photutils.centroids.gaussian]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg\n",
      "    (277.73276576, -10.49944634)>, <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg\n",
      "    (277.73153159, -10.49993895)>, <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg\n",
      "    (277.73176784, -10.49978128)>]\n",
      "277.7320217328176\n",
      "0.0005348774944543936\n",
      "Average RA: 277.7320217328176 ± 0.0005348774944543936\n",
      "Average Dec: -10.499722188279634 ± 0.00020540453144858676\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Input data contains non-finite values (e.g., NaN or inf) that were automatically masked. [photutils.centroids.gaussian]\n"
     ]
    }
   ],
   "source": [
    "# Measuring the centroid position for each map\n",
    "Positions684 = []\n",
    "Positions684_values = np.empty((len(ListCrops684),2))\n",
    "for i in range (0, len(ListCrops684)):\n",
    "    Position = photutils.centroids.centroid_2dg(ListCrops684[i].data, error=None, mask=None)\n",
    "    w = wcs.WCS(ListCrops684[i].header)\n",
    "    PositionRADec = w.pixel_to_world(Position[0],Position[1])\n",
    "    Positions684.append(PositionRADec)\n",
    "    # RA values\n",
    "    Positions684_values[i,0] = PositionRADec.ra.deg\n",
    "    # Dec values\n",
    "    Positions684_values[i,1] = PositionRADec.dec.deg\n",
    "\n",
    "print(Positions684)\n",
    "\n",
    "# Find the average and standard deviation in RA and Dec\n",
    "print(np.average(Positions684_values[:,0]))\n",
    "print(np.std(Positions684_values[:,0]))\n",
    "\n",
    "# Find the average and standard deviation in RA and Dec\n",
    "RA684_av = np.average(Positions684_values[:,0])\n",
    "RA684_std = np.std(Positions684_values[:,0])\n",
    "Dec684_av = np.average(Positions684_values[:,1])\n",
    "Dec684_std = np.std(Positions684_values[:,1])\n",
    "\n",
    "print('Average RA: '+str(RA684_av)+' ± '+str(RA684_std))\n",
    "print('Average Dec: '+str(Dec684_av)+' ± '+str(Dec684_std))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Input data contains non-finite values (e.g., NaN or inf) that were automatically masked. [photutils.centroids.gaussian]\n",
      "WARNING: Input data contains non-finite values (e.g., NaN or inf) that were automatically masked. [photutils.centroids.gaussian]\n",
      "WARNING: Input data contains non-finite values (e.g., NaN or inf) that were automatically masked. [photutils.centroids.gaussian]\n",
      "WARNING: Input data contains non-finite values (e.g., NaN or inf) that were automatically masked. [photutils.centroids.gaussian]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg\n",
      "    (277.68266508, -10.3689247)>, <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg\n",
      "    (277.68269599, -10.36885723)>, <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg\n",
      "    (277.68261358, -10.36883983)>, <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg\n",
      "    (277.68252566, -10.3689537)>]\n",
      "Average RA: 277.6826250783497 ± 6.450551699358738e-05\n",
      "Average Dec: -10.368893864282477 ± 4.6887071628262587e-05\n"
     ]
    }
   ],
   "source": [
    "# Measuring the centroid position for each map\n",
    "Positions685 = []\n",
    "Positions685_values = np.empty((len(ListCrops685),2))\n",
    "for i in range (0, len(ListCrops685)):\n",
    "    Position = photutils.centroids.centroid_2dg(ListCrops685[i].data, error=None, mask=None)\n",
    "    w = wcs.WCS(ListCrops685[i].header)\n",
    "    PositionRADec = w.pixel_to_world(Position[0],Position[1])\n",
    "    Positions685.append(PositionRADec)\n",
    "    # RA values\n",
    "    Positions685_values[i,0] = PositionRADec.ra.deg\n",
    "    # Dec values\n",
    "    Positions685_values[i,1] = PositionRADec.dec.deg\n",
    "\n",
    "print(Positions685)\n",
    "\n",
    "# Find the average and standard deviation in RA and Dec\n",
    "RA685_av = np.average(Positions685_values[:,0])\n",
    "RA685_std = np.std(Positions685_values[:,0])\n",
    "Dec685_av = np.average(Positions685_values[:,1])\n",
    "Dec685_std = np.std(Positions685_values[:,1])\n",
    "\n",
    "print('Average RA: '+str(RA685_av)+' ± '+str(RA685_std))\n",
    "print('Average Dec: '+str(Dec685_av)+' ± '+str(Dec685_std))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative pixel error for Flight 684 in RA: 0.4231997758320439\n",
      "Relative pixel error for Flight 684 in Dec: 0.16251787103624302\n",
      "---\n",
      "Relative pixel error for Flight 685 in RA: 0.051037332126793955\n",
      "Relative pixel error for Flight 685 in Dec: 0.03709746326631732\n"
     ]
    }
   ],
   "source": [
    "# Pixel scale in Band E\n",
    "BandE_pix = 0.0012638888888889\n",
    "print('Relative pixel error for Flight 684 in RA: '+str(RA684_std/BandE_pix))\n",
    "print('Relative pixel error for Flight 684 in Dec: '+str(Dec684_std/BandE_pix))\n",
    "print('---')\n",
    "print('Relative pixel error for Flight 685 in RA: '+str(RA685_std/BandE_pix))\n",
    "print('Relative pixel error for Flight 685 in Dec: '+str(Dec685_std/BandE_pix))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Difference in RA: -0.049396654467898315\n",
      "Difference in Dec: 0.13082832399715727\n",
      "These differences needs to be added to the CRVAL values in each extension of the Level 2 data headers\n"
     ]
    }
   ],
   "source": [
    "# Measuring the difference\n",
    "DeltaRA = RA685_av - RA684_av\n",
    "DeltaDec = Dec685_av - Dec684_av\n",
    "\n",
    "print('Difference in RA: '+str(DeltaRA))\n",
    "print('Difference in Dec: '+str(DeltaDec))\n",
    "print('These differences need to be added to the CRVAL values in each extension of the Level 2 data headers')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PoLiteWIP",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
